Persistence exponents in a 3D symmetric binary fluid mixture 
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The persistence exponent, 9, is defined by Nf ~ t , where t is the time since the start of the 
coarsening process and the "no-flip fraction", Nf, is the number of points that have not seen a 
change of "color" since t = 0. Here we investigate numerically the persistence exponent for a binary 
fluid system where the coarsening is dominated by hydrodynamic transport. We find that Nf 
follows a power law decay (as opposed to exponential) with the value of 9 somewhat dependent on 
the domain growth rate (L ~ t" , where L is the average domain size), in the range 9 = 1.23 ± 0.1 
(q = 2/3) to 9 = 1.37 ± 0.2 (q = 1). These a values correspond to the inertial and viscous 
hydrodynamic regimes respectively. 

PACS numbers: 64.60.Ht, 64.75+g, 82.20.Wt 



I. INTRODUCTION 

Persistence exponents can be defined for systems where 
there is an order parameter whose time evolution can be 
followed at each point, x, in the system. It is most easily 
understood for an order parameter, 0(x) that takes just 
two equilibrium values (spin up/down), but the concept 
is easily generalised (e.g. red/blue fluid for binary fluid 
mixtures). The density of points for which sgn[(/)(x) — ((/>)] 
has not changed sign up to time t as the system coarsens, 
the "no-flip fraction", Np, will, for a conserved order pa- 
rameter, or a quench to zero temperature, typically decay 
as a power law, Np ~ t~^ , where 9 is the persistence ex- 
ponent. 

Persistence exponents were first investigated by Bray, 
Derrida and Godreche ||l|,H, in the context of one di- 
mensional diffusive systems. For example, for the g-state 
Potts model in one dimension, Derrida et. al. B found 
from simulations, and later proved analytically [0 , that 
e = 3/8 for q = 2, 61 ~ 0.53 for q = 3 and 6* ^ 1 as 
q — > oo. It is also possible to define persistence expo- 
nents for systems with non-conserved order parameters 
at finite temperature by coarse-graining M. 

As with most critical systems, analytical calculations 
are difficult; a few exact results exist [^|,^. Some mean 
field calculations have been done [^, however, mean field 
approximations also predict that 9 is not independent of 
the other critical exponents. Mean field theory assumes 
that the order parameter dynamics is a Markov process, 
and Derrida et. al. and Majumdar et. al. iGUq] argue 
that in general this is not the case, and thus 9 is, in fact, 
a new critical exponent independent of the two static and 
two dynamic exponents already known. 

In this paper we present an investigation of the persis- 
tence behavior of our simulation of 3D spinodal decom- 
position in a binary fluid mixture, other results of which 
are reported elsewhere [QJ^] . Brief details of the theoreti- 
cal model and simulation method are given first, followed 



in Section III by a theoretical approach to persistence be- 
havior in this system. In Section IV the numerical anal- 



ysis is described, with the results for the flip rate given 
in Section Mand estimates of the persistence exponent 
in Section W\ Finally, in Section VII , the results are 
summarised. 



II. MODEL AND SIMULATION 

A symmetric binary fluid mixture differs from a purely 
diffusive system such as an alloy in that the phase sepa- 
ration is assisted by hydrodynamics. There is an initial 
diffusive period during which an interlocking structure of 
single-fluid regions is formed, separated by sharp inter- 
faces. The interfaces then take over as the driving force, 
displacing the fluid as they flatten and shrink, leading 
to a much more rapid domain growth rate of L ~ i (in 
3D) compared to diffusive growth of L ~ t^^^, where L is 
the average domain size. This linear growth was first pre- 
dicted by Siggia 19] . Furukawa |lfl| later predicted that as 
the inertial effects came into play, the growth rate would 
slow to i ^ t^/^ . These growth rates have been observed 
in simulation [0 11 ] and linear growth has been observed 
experimentally 12]. Recent suggestions that the growth 
rate may slow still further to L ^ t^/'^ (or slower) ]|l^ are 
not supported by our simulation work ]R] or theory WM , 
and will not concern us here. 

For numerical work, we use the following model free 
energy: 

F = jdv\^-^4>' + f 0' + /51np+ |]V0]2| , (1) 

in which A, B and k are parameters that determine the 
interfacial width (^ — y/K/2A), and interfacial tension 
(ct — y/8KAy9B^); </> is the usual order parameter (the 
normalized difference in number density of the two fluid 



species); p is the total fluid den sity , which remains (vir- 
tually) constant throughout |ll^,|l6| . We choose A/ B = 1 
so (/) = ±1 in equilibrium. Our simulation [ [L7[ uses a 
lattice Boltzmann (LB) method llSlJlSfl with a cubic lat- 
tice with nearest and next-nearest neighbor interactions 
(D3Q15). It was run on Cray T3D and Hitachi SR-2201 
parallel machines with system sizes up to 256^. 

The data used for the analysis of persistence behav- 
ior corresponds to that used in our earlier studies of the 
growth exponents ||7|] . Details of the simulation parame- 
ters for each run, and selection of usable data can also be 
found in 0]. The filters used to eliminate diffusive and 
finite size effects mean that the good data from any single 
256'^ run lies within 20 < L < 64 in units of the lattice 
spacing. The data was fitted to L = h{t — iint)", where b 
is a prefactor determined by the physical parameters and 
tint is a non-universal adjustment to the zero point on the 
time scale dependent on the initial diffusive period. Fig- 
ure H shows two sample runs and the fitted curves. The 
values of a and iint will be used in the subsequent anal- 
ysis of the persistence behavior. By using characteristic 
length (Lq = r]^ /{pa)) and time (ip = V^/{p<^'^)) scales 
uniquely defined by the physical parameters, the data 
from all the runs was combined into a single L, t plot 
covering a linear region, through a broad crossover, to 
t^/"^ 0- The value of the growth exponent, a, where the 
domain size, L ~ i", was found to range from 1.0 to 0.67, 
thus the data spans the full range from viscous to inertial 
hydrodynamic growth. The breadth of the crossover re- 
gion justifies the use of a single effective exponent, a, to 
fit any single run. It is convenient to use the value of Lq 
(in simulation units) to distinguish the different simula- 
tion runs. The viscous regime corresponds to Lq > 1) the 
inertial regime to Lq < 0.001, with intermediate values 
corresponding to the crossover region. 



III. FLIP-RATE MODEL 

For the 3D binary fiuid system it is unrealistic to ex- 
pect to derive anything but the simplest approximate 
results using a theoretical approach. However, this will 
serve to illuminate the persistence-related quantities un- 
der discussion. The simple results to be derived in this 
section will then be of assistance in the numerical analy- 
sis that follows. 

Consider the step by step {At = 1 in simulation units) 
time evolution of the value of the order parameter at a 
single point as the system phase separates and coarsens 
from the initially completely mixed state. A typical ex- 
ample is represented schematically in Fig. H^. In order 
to focus on the behavior in the hydrodynamic regime, 
an initial state is chosen at a reference point, istart, a 
time corresponding to when domains have coarsened to 
a size L ^ Lmin that marks the onset of purely hydrody- 
namic behavior (determined by the point at which diffu- 
sive growth has fallen below 2% R). The "no-fiip frac- 



tion", Np{t /tstart) is then defined as the fraction of sites 
that have not changed color since tstart- Scaling by igtart 
is equivalent to choosing units such that igtart — 1, which 
is in any case the value of the initial reference time used 
in other studies. 

In order to derive an approximate functional form for 
Np{t/tstart), it is useful to define two further quantities, 

• the flip rate, PF{t) is simply the proportion of sites 
that changed color between time step t — 1 and t, 
and, 

• the flip probability, P{t,ti), is the probability that 
a site changes color at time i, given that it last 
changed color at time ti. 

These two quantities are related as follows: 



Pp{t) 



dhPF{h)P{t,ti] 



(2) 



where the sum over discrete time in our simulations has 
been approximated by an integral. Equation (0) says 
that the flip rate at time t is given by all the points that 
last flipped at time ii and are due to flip again at time t, 
i.e. P{t,ti), integrated over all possible prior flip times, 
tint < ti < t, and weighted by the number of sites with 
prior flip time ti, i.e. Pp{ti). The lower limit of the 
integral is set to tint because, as can be seen from Fig. 
0, tint corresponds to the natural zero point on the time 
scale, the time at which, (ignoring diffusion and the flnite 
width of the interfaces) , the domain size would be zero so 
that, in effect, every site flipped at t = tint- Henceforth 
tint will be usually be set to zero, to simplify the algebra. 

Equation (0) is an integral equation in two unknown 
functions, P(t,ti) and PF{t). A solution for P(t,ti) 
can be obtained by making some assumptions about the 
asymptotic form of the simpler quantity, Pp{t), for the 
3D hydrodynamic system. The dynamics are determined 
by the basic scaling growth law, L — b{t — tint)", where 
the prefactor v depends on the system parameters (den- 
sity, viscosity, interfacial tension). Once the initial dif- 
fusive period is over, and the domains of red and blue 
fluid are separated by well-formed interfaces, lattice sites 
change color from red to blue when an interface moves 
across them as the domains grow in size. The flip rate, 
PF{t), is thus given by the rate at which the interface 
moves through the system sweeping out a volume of 
points that change color. So, an estimate of the area 
of the interface, and its average speed, will yield an esti- 
mate for Pf (t) . This is the "flip-rate model" . 

The area of the interface is given approximately by 
A{t) = clV/L where V is the system volume and cl is 
a prefactor of order unity. As the system coarsens, the 
interfaces must move to accommodate the enlargement 
of the domains. It has been shown B that the average 
fiuid velocity is comparable in magnitude with dL/dt, 
so it is reasonable to estimate the speed of the interface 
as Cy di/dt, where Cy is another prefactor of order unity. 



The volume swept out per unit time will be A{t) Cy dL/dt. 
Combining both cl and Cy into a single prefactor, c, the 
flip rate per unit volume will be approximately, 



Pp{t) 



1 dL 
'L~dt 



ca 



(3) 



This approximation for Ppit) diverges at t = 0. In the 
discrete time of our simulations, we start with timestep 
t = 1 so no divergences arise; in the continuum approxi- 
mations it is appropriate to start at time i = 0, thus we 
will need to take care that no integrals diverge at their 
lower limits. 

Substituting Ppit) — ca/t into Eq. (0) gives. 



/■* t 

/ dii -P(t, ii) = l, 

Jo *i 



where have set ijnt = 0. Equation 
P{t,ti) of the form, 

P(t,ti) = (/3-l)tf-V^ 



(4) 
has solutions for 



/3>1, 



(5) 



where /3 is an arbitrary exponent, as can readily be veri- 
fied by substitution. We have P{t,ti) ^ for ii ^ so 
the integral is well-behaved at ii = 0. The condition. 



/■oo 

/ dtp{t,h)^l, 

Jtl 



(6) 



{i.e. all sites do, eventually, flip) is also satisfled, so 
P{t,ti) is a properly normalised probability. 

An expression can now be written down for the "no 
flip fraction", Npit/tstait), which is the fraction of sites 
that have not changed color since istart, 



Np = 



i start 



dh / dt2P{t2,ti)PF{tl), 



(7) 



i.e. we count every point whose last flip was before istart 
and whose next flip is after time t. Substituting for Pp 
and P{t2,ti) gives. 



tsta 



Nf ^ca I dii (/3 - 1)*^' 
/o 

e 



P-2 



dt2t^^, 



ca I t 



'-start 



t 



61 > 0, 



(8) 



setting 61 = (/?— 1) as the persistence exponent. The flip- 
rate model has thus provided an expression for Np that 
has the exponent, 0, in the prefactor as well as being 
the asymptotic power of the decay. This puts helpful 
constraints on the data analysis, although for it to be 
really useful, the prefactor, c, in Eq. (||) has to be pinned 
down as precisely as possible. In a log- log plot. 



logNp = log(ca/6l) - 6'log(i/istart). 



(9) 



so plotting Nf against i/istart should give a straight line 
when viewed log-log, with slope ~9, for t ^ istart- The 



initial data near t > tstart will depart from this asymp- 
totic behavior and tend to Np = 1 because by definition 
nothing flipped yet for t ~ istart- The fit to the asymp- 
totic slope should, however, have an intercept at ca/6 for 

^ — ^start • 

It is possible to go on to define a hierarchy of persis- 
tence quantities, of which Np is the first Q. The "one flip 
fraction" , Op is the density of points that have changed 
color exactly once since time tstart- Within the flip-rate 
model, an expression for Op for this system can be de- 
rived in the same way as for Np. 

Of (t/istart) = 

oo /■*2 /"^start 

dt2 / dtf / dhP{t2,tf)P{tfM)PF{tl). 

t JtBtart Jo 

(10) 

The "one flip" occurs at time tf with tgtart < i/ < t2- 
The probability of an initial flip at ti with < ti < tgtart, 
then flipping at time tf is P{tf,ti). The probability of 
flipping again at time t2 is then P{t2,tf)P{tf,ti), pro- 
vided the two flip events are independent of each other 
(so the probabilities can be multiplied). In other words, 
this model assumes that the flip probability, P{t2,tf) is 
only dependent on the last time the site flipped, it has no 
memory of earlier flips. This requires a Markov process 
for the coarsening dynamics, which Derrida et. al. and 
Majumdar et. al. |^,g claim is not generally the case 
for the non-hydrodynamic systems they have considered. 
This therefore constitutes another untested assumption 
in the flip-rate model described here. 

Substituting for Pp and P{t2,ti) as before gives. 



2 / ^i-f 



Op = / di 

COi ( istart 
t 



dti 



^sta 



COL „ ^t 

dt^ — e ^ 



tl t«+l 



9+1 



1-ln 



t 



(11) 



after integrating by parts. This has a dominant loga- 
rithmic term so doesn't have an asymptotic power law 
decay. (In contrast, in the ID Ising model, Op appears 
to have the same asymptotic behavior as Np, see P,pi.) 
This could, of course, be a result of the limitations of the 
flip-rate model, however, the simulation results suggest 
that the asymptotic behavior of Op is different from that 
of Np (see Figs, y and below). Again, this expression 
is only valid for t ^ istart since Op — > for t — > tgtart- 
As with Np, the (same) critical exponent Q appears both 
as an exponent and as a prefactor so providing an extra 
constraint on any flts to simulation data. 

Nothing in the theory so far is specific to three space 
dimensions. However, in two dimensional binary fiuid 
systems, the interface doesn't completely interconnect 
|19| , so the approximation Pp = cajt is unlikely to work 
without modification. The two dimensional form of the 
theory will not be pursued further here. 



IV. NUMERICAL ANALYSIS 



V. FLIP RATE 



The flip-rate model just described will now be used to 
guide the analysis of the data from our spinodal decom- 
position simulations. The data from the largest (256'^) 
runs is typically saved every 200 to 500 time steps, 
Aistrido- It is also coarse-grained from a 256'^ lattice 
down to 128'^ by averaging over groups of eight neigh- 
boring lattice points. This data can be used to calculate 
the four key quantities discussed in the previous section, 
Pp, P{ti,t2), Np and Op, provided the spatial and tem- 
poral coarse-graining makes no significant difference to 
the results. 

To investigate temporal coarse-graining, a special run 
on a 96^ grid was done with the key quantities calculated 
every time step. This extra computation significantly 
slows down the computer run time so it isn't practical to 
perform larger simulations with single-step calculations. 
The results from a single-step calculation were then com- 
pared with the same run analysed over a stride of 100 
time steps, see Fig. 0. The result of this comparison 
showed excellent agreement between the single-step and 
strided estimates of Pp (t) . 

In fact, the errors introduced by analysing only ev- 
ery Atstridc steps rather than every single step arise from 
miss-counting sites that flip more than once in that pe- 
riod. If the flip rate, Ppit), is small then the number of 
multiple flips will be much smaller still so the errors will 
negligible. In all runs analysed, Pp, was found to be small 
even over the largest Aigtrido used, i.e. Pp x Aistridc ^ 1 
(a value of one would mean every site flipped). 

The shape oi Pp{t) in Fig. 0for times t < tstart clearly 
shows a transition from diffusive to hydrodynamic be- 
havior in the initial fall (diffusion), rise (interfaces start 
to move) and fall again (hydrodynamic). Thus Pp{t) 
is a sensitive indicator of the system dynamics, and the 
choice of tstart = 1500 in this %& system is confirmed to 
be located, as desired, near the beginning of the hydro- 
dynamic regime but not too close to the transition from 
diffusive behavior. 

The effect of spatial coarse-graining was checked by 
comparing the results of analysing strided data from 128^ 
and 256^ runs with identical parameters apart from sys- 
tem size. Excellent agreement was found everywhere ex- 
cept in the tail of P(i, tgtart), which can be explained by 
poor statistics in this region. The conclusion, therefore, 
is that spatial and temporal coarse-graining on the scale 
used in this analysis does not introduce significant errors 
into the results. 

From the main simulation data, good statistics for Pp, 
Np and Op can be obtained, as every lattice point can 
be assessed according to if/ when it next fiips after tstart- 
The statistics are not so good for estimating P(i,tstart), 
because the reference base is reduced to only those sites 
that flip between igtart and Aigtrido rather than the whole 
system; recall that P(i,istart) refers to the first flip hap- 
pening exactly at igtait- 



In order to make full use of the model presented in 
Section [II , it is necessary to investigate whether the flip 
rate really follows the approximate theoretical expres- 
sion, Pp — ca/t, Eq. (0), and whether the prefactor, c, 
can be evaluated sufficiently accurately from the simula- 
tion data to allow a strongly constrained fit to be done 
for 6 in the analysis of Np and Op . Figure ^ shows Pp 
for each simulation run in a log- log plot, with \/t and 
1.5/i also shown for comparison. The prefactor, c, clearly 
varies somewhat over time, and with a, although the 
variation with time can partly be explained since early 
times, t ^ tstart, are expected to differ from asymptotic 
behavior. 

To determine the values and variation in c more ac- 
curately, a linear plot of the flip rate, Pp, scaled by 
{t — iint)/a is shown in Fig. 1^. It can be seen that the 
results roughly split into two groups corresponding to the 
runs found to be in the viscous and lower crossover re- 
gions (Lq ^ 0.1), with prefactors around c ~ 1.25, and 
runs in the inertial regime {Lq < 0.01), with c ~ 2.2. 

If this difference in the value of c is real, it implies 
that the geometry of the separating domains is different 
between the viscous and inertial regimes. Any such dif- 
ference is not apparent to the eye in visualisations of the 
interface pO]. There is some evidence from analysis of 
the structure factor Q that there is a structural differ- 
ence in the domain structure, but it is unclear whether 
this is enough to account for the observed difference in 
the value of c. However, there is another possible inter- 
pretation of the data, which is that the flip rate actually 
exceeds ca/t for the inertial regime runs. For, once the 
viscosity is low enough for inertial effects to become sig- 
nificant, the interface may exhibit capillary excitations. 
This oscillation of the interface could increase the mea- 
sured flip rate as sites near the interface repeatedly flip 
back and forth, but without contributing to the sweep- 
ing that removes sites from the "no flip" category. If 
capillary excitations have a significant effect on the ffip 
rate, this part of the flip rate needs to be discounted in 
the subsequent fitting to determine the persistence ex- 
ponent. Since the extent of any excess flip rate cannot 
be determined at this stage, a variable prefactor, c, has 
therefore been carried through the fitting process for 9. 



VI. PERSISTENCE EXPONENT 

The "no fiip fraction", Np, is the primary quantity 
of interest in this work, being the quantity from which 
the persistence exponent 9, is defined, Np ~ t~^ . The 
simulation results for Np, with the time scaled as {t — 
^int)/istart and points for times at which L < Z/min or 
L > 64 removed, are plotted in Fig. g. The runs in the 
inertial regime have tails that fall to around Np = 0.08, 
while the viscous runs stop (reach domain size L = 64) 



TABLE I. Summary of results for c and 9. 



Lo 


a 


c 




e 


0.0003 


0.67 


1.5- 


1.9 


1.20- 


- 1.23 


0.00095 


0.67 


1.5- 


2.2 


1.20- 


- 1.32 


0.01 


0.75 


1.3- 


1.9 


1.25- 


- 1.40 


0.15 


0.80 


1.0- 


1.5 


1.15- 


- 1.24 


0.95 


0.95 


0.85 - 


- 1.4 


1.20- 


- 1.40 


5.9 


1.0 


0.85 - 


- 1.4 


1.25- 


- 1.50 


36 


1.0 


0.85- 


- 1.3 


1.30- 


-1.55 



at around Np — 0.15. There is thus barely one decade 
within which to estimate the slope of the tail. Although 
a simple linear fit cannot produce a very accurate result, 
superficial inspection suggests that is slightly greater 
than unity, and fairly reproducible. 

If we could use a fixed value of the prefactor, c, with 
confidence, the fit could be constrained to intercept at 
ca/9, Eq. (||). Instead, a family of fits was done P0[ , 
covering a range of values of the prefactor, c, and these 
were compared with a similar family of fits for the one 
flip fraction. Op- As can be seen in Fig. ^, Op at best 
falls to only 0.2 by the end of the simulation data, so can- 
not be regarded as having reached asymptotic behaviour 
yet. The Op data was therefore used only to confirm the 
direction in which Np is tending towards its asymptote. 
From this fitting procedure, a range of consistent values 
of both c and 9 were obtained, these are summarised in 
Table |. 

A sample plot for the run with Lq ~ 0.0003 with Np 
and the fitted expression in Eq. (|^), with high, median 
and low values of c and 6, is shown in Fig. H. 

Fits to P{t, istart) = 0/t{tsti,n/ty wcrc also done, pro- 
ducing values of 9 consistent with those determined from 
Np and Op. Since P(t, tstart) does not depend on c or a, 
this is a useful check for consistency, but because of the 
poorer statistics for P(i, istart), no improvement in the 
range of values for 9 was obtained by this extra step. 

The complete set of estimates for 9 from Table H is 
shown against the growth exponent, a in Fig. H with 
error bars indicating the range of values obtained The 
results are consistent with the persistence exponent, 9, 
lying somewhere around 1.37 for the viscous regime and 
1.23 for the inertial regime. A single value oi9 = 1.333 is 
also just about consistent with the data, although some 
variation with a seems more likely. 

The prefactor, c, could be as low as 1.0 for viscous 
runs and rises to nearer 2.0 for inertial regime. How- 
ever, a value of c — 2.2 as implied by the measured fiip 
rate for the inertial regime appears to be inconsistent 
with the persistence data for Np, Op, and P{t, tstart)- 
The flip rate is thus sufficiently sensitive to distinguish 
between viscous and inertial regime dynamics, and sug- 
gests the presence of some mechanism such as capillary 
waves in the inertial regime that raises the measured flip 
rate above that predicted by the growth of the domains. 



VII. CONCLUSIONS 

In this investigation of the persistence exponents for a 
3D hydrodynamic spinodal system, it has been observed 
that the decay of the "no-flip fraction", Np, is a power 
law (as opposed to exponential). It has not been pos- 
sible to produce a very accurate determination of the 
persistence exponent, 9, but limits have been placed on 
the likely value, and it has been shown that some varia- 
tion with the effective growth exponent, a, as one crosses 
slowly from viscous (a = 1) to inertial (a — 2/3) coars- 
ening, should be allowed for in any future studies. The 
best estimate of 9 for the viscous regime {a = 1), is 
9 = 1.37 ± 0.2. From the point of view of fundamental 
critical exponents, the value of 9 in the inertial regime is 
perhaps more significant, since this is (at least on current 
evidence) the asymptotic behavior for spinodal decompo- 
sition, and here the best estimate is 9 = 1.23±0.1. These 
represent the first estimates of in a system dominated 
by hydrodynamics rather than diffusion. However, the 
error limits quoted here do not allow for systematic er- 
rors arising from the use of the flip rate model itself. More 
accurate determination from simulations would likely re- 
quire significantly larger system sizes, which is difficult 
to envisage at present computing power. 
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FIG. 2. L vs. t for the runs with Lq — 5.9 (circles) and 0.0003 (diamonds). The region used for fitting is delimited by 
Lmin < L < Lmax ~ 64, aud the fits (solid) (to a = 1 and a — 2/3 respectively) are projected back to show the intercepts, tint- 
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FIG. 3. Flip rate, Pf, for 96"' system single step (dots) and stride=100 (crosses) 
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FIG. 4. Flip rate, Pf, for simulation runs as identified by Lo values in key, log-log plot. Also shown for comparison, 1/t and 
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FIG. 6. No-flip fraction, Nf, for all runs. 
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FIG. 8. Fitted curves for Nf for the run with Lq = 0.0003. The no flip fraction, Nf, (squares) calculated from the numerical 
data is shown with the upper, median and lower fits to Eq. (H) using c — 1.9, 6 — 1.23 (dashed) c — 1.7, 9 — 1.215 (solid) and 
c — 1.5, 6 = 1.21 (dotted) respectively from Table |. 
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